1. Coping with Mortality

Learning Objectives

Understand options for modeling background mortality based on life tables and parametric mortality models.
Construct a cause-deleted life table to model cause-specific and non-cause-specific death.

Approaches to Modeling Mortality

Standard approach is to draw from life table data:

age	lx	dx	mx	qx
0	100000	553.22	0.00556	0.00553
1	99447	38.18	0.00038	0.00038
2	99409	23.16	0.00023	0.00023

lx is the total number living in a hypothetical cohort (of 100,000).
dx is the number of deaths in the age interval.

Approaches to Modeling Mortality

Standard approach is to draw from life table data:

age	lx	dx	mx	qx
0	100000	553.22	0.00556	0.00553
1	99447	38.18	0.00038	0.00038
2	99409	23.16	0.00023	0.00023

mx is the death rate within the interval.
qx is the probability of death within an interval.

Life Table Data

Death probabilities are useful for discrete time Markov models.

age	lx	dx	mx	qx
0	100000	553.22	0.00556	0.00553
1	99447	38.18	0.00038	0.00038
2	99409	23.16	0.00023	0.00023

Life Table Data

Death rates are useful for discrete event simulation (DES) models.

age	lx	dx	mx	qx
0	100000	553.22	0.00556	0.00553
1	99447	38.18	0.00038	0.00038
2	99409	23.16	0.00023	0.00023

Life Table-Based Mortality

A common approach is to use a life table-based lookup table in the model:

tpDn_lookup <-
    c("(34,44]" = 0.0017,  
      "(44,54]" = 0.0044,
      "(54,64]" = 0.0138,
      "(64,74]" = 0.0379,
      "(74,84]" = 0.0912,
      "(84,100]" = 0.1958)

Highlighted row: annual probability of death for 35-44 year old.

Life Table-Based Mortality

What if age bins are coarse (e.g., 0, 5, 10,… vs. 0,1,2,3,…)?
- Assume death probability/rate is constant within age bins (and cycles)?
- In a DES model, we need to sample a time to death.
Can sample from survival function in underlying life table.
Ages in life table often topcoded (e.g., 85, 100, etc.)

Alternative: Parametric Mortality Model

Fit a parametric mortality model to the life table data.
Reduces dimensionality to a few parameters.
Use these parameters to sample time of death, obtain mortality rate/probability at any arbitrary time, etc.

Obtaining Life Table Data

The R packages MortalityLaws and demography facilitate life table data aquisition and modeling.

Common sources:

Processing Life Table Data

Code below extracts U.S. life table data from the human mortality database (HMD).

Processing Life Table Data

Code below extracts U.S. life table data from the human mortality database (HMD).

hmd_usa <-
    demography::hmd.mx("USA",
                       username = "<<your user name>>",
                       password = "<<your password>>", "USA")

Processing Life Table Data

Code below processes HMD data and constructs a (2019) life table from it.
Radix (cohort size) is 100,000 individuals.

Processing Life Table Data

Code below processes HMD data and constructs a (2019) life table from it.
Radix (cohort size) is 100,000 individuals.

mortality_year = 2019 
radix = 100000

lt = 
    hmd_usa %>% 
    demography::lifetable(.,series = "total", years = mortality_year) %>% 
    as_tibble() %>% 
    rename(age = x) %>% 
    mutate_at(vars(lx,dx), function(x) x * radix)

US Life Table (2019)

Highlighted column: Remaining life expectancy at exact age x.

age	lx	dx	ex	qx	mx
0	100000	553.221	78.963	0.00553	0.00556
1	99447	38.180	78.402	0.00038	0.00038
2	99409	23.160	77.432	0.00023	0.00023
3	99385	17.689	76.450	0.00018	0.00018

Alive-Dead Model

As a refresher, let’s construct a very simple alive-dead model using life table death probability (qx) values.

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

No specific strategies or policies under consideration; we’re just trying to replicate mortality using life table data for a hypothetical cohort of 100,000 infants.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )

Discrete time will be in annual cycles.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )

The lifetable “strategy” will draw on the HMD life table data processed earlier.

Age-Dependent Transition Matrix

We’ll next construct a transition probability matrix using the age-specific probability of death field (qx) from the life table data.
To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.

Age-Dependent Transition Matrix

Life Table:

lt[1,c("age","qx")]

# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1     0 0.00553

lt[90,c("age","qx")]

# A tibble: 1 × 2
    age    qx
  <dbl> <dbl>
1    89 0.123

Transition Matrix:

params$mP[[1]][,,"lifetable"]

       to
from    Alive   Dead     
  Alive 0.99447 0.0055322
  Dead  0       1

params$mP[[90]][,,"lifetable"]

       to
from    Alive   Dead   
  Alive 0.87677 0.12323
  Dead  0       1

Construct the Markov Trace

sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}

Construct the Markov Trace

Output: named list object with full Markov trace for each state.

Construct the Markov Trace

Output: named list object with full Markov trace for each strategy.

trace <- 
  sim_cohort(params)
trace[["lifetable"]][1:10,]

   Alive   Dead
0 100000   0.00
1  99447 553.22
2  99409 591.40
3  99385 614.56
4  99368 632.25
5  99354 646.46
6  99340 659.67
7  99328 671.69
8  99317 682.52
9  99307 693.14

Checkpoint: Life Expectancy at Birth

Life Expectancy

We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
According to the underlying life table data, life expectancy at birth is 78.96315 years.
Let’s verify we can replicate this with our Markov model.

Life Expectancy

We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
According to the underlying life table data, life expectancy at birth is 78.96315 years.
Let’s verify we can replicate this with our Markov model.

age	qx	ex
0	0.00553	78.963

Life Expectancy

Define the payoff for life expectancy (1 if alive, 0 otherwise)

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

Calculate life years in each cycle.¹

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

Define a function for cycle adjustment.¹

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

Sum of values (with cycle adjustment applied) is total life expectancy.

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

Age	Life Expectancy (Life Table data)	Life Expectancy (Markov-Life Table)
0	78.963	78.924

We can repeat this exercise for different cohort starting ages and verify things line up.

Life Expectancy

Mortality Modeling

Depending on context and modeling environment, you may need to model mortality.
The MortalityLaws package has a number of mortality models we can draw from:

Mortality Modeling

Name	Model	Code
Infant mortality
Opperman	$\mu[x] = A/\sqrt{x} - B + C\sqrt{x}$	opperman
Weibull	$\mu[x] = 1/\sigma * (x/M)^{(M/\sigma - 1)}$	weibull
Accident hump
Inverse-Gompertz	$\mu[x] = [1-\exp((M-x)/\sigma)]/[\exp((M-x)/\sigma)-1]$	invgompertz
Inverse-Weibull	$\mu[x] = 1/\sigma * (x/M)^{[-M/\sigma - 1]} / [\exp((x/M)^{(-M/\sigma)}) - 1]$	invweibull

Name	Model	Code
Adult mortality
Gompertz	$\mu[x] = A \exp[Bx]$	gompertz
Gompertz	$\mu[x] = 1/\sigma \exp[(x-M)/\sigma]$	gompertz0
Makeham	$\mu[x] = A \exp[Bx] + C$	makeham
Makeham	$\mu[x] = 1/\sigma \exp[(x-M)/\sigma)] + C$	makeham0
Perks	$\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]$	perks
Strehler-Mildvan	$\mu[x] = K \exp[-V_0 (1-Bx)/D]$	strehler_mildvan
Adult and/or old-age mortality
Van der Maen	$\mu[x] = A+Bx+Cx^2 + I/[N-x]$	vandermaen
Beard	$\mu[x] = \exp(Bx)/[1+KA \exp(Bx)]$	beard
Beard-Makeham	$\mu[x] = A \exp(Bx)/[1 + KA \exp(Bx)]$	beard_makeham
Gamma-Gompertz	$\mu[x] = A \exp(Bx)/(1+AG/B[\exp(Bx)-1])$	ggompertz

Name	Model	Code
Old-age mortality
Van der Maen	$\mu[x] = A + Bx + I/[N-x]$	vandermaen2
Quadratic	$= A + Bx + Cx^2%	quadratic
Kannisto	$\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]$	kannisto
Kannisto-Makeham	$\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]+C$	kannisto_makeham

Name	Model	Code
Full age range
Thiele	$\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}$	thiele
Wittstein	$q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}$	wittstein
Heligman-Pollard	$q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x$	HP
Heligman-Pollard	$q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]$	HP2
Rogers-Planck	$q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}$	rogersplanck
Martinelle	$\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}$	martinelle
Carriere	$l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)$	carriere1
Kostaki	$q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x$	kostaki

Mortality Modeling

Another nice feature of the package is that each mortality model type has its own defined function to calculate the mortality rate (or probability) given an age and model coefficients.

Mortality Modeling

Generally speaking, we need three inputs:

age: Ages for lifetable
dx: The number of deaths between exact ages x and x+1.
lx: Number of survivors to exact age x.

Gompertz Model

Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
Gompertz reduces mortality to two parameters.

Gompertz Model

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

gom_fit <- MortalityLaw(
    x  = ages,      # vector with ages
    Dx  = deaths,   # vector with death counts
    Ex  = exposure, # vector containing exposures
    law = "gompertz",
    opt.method = "LF2")

Gompertz Model

# Fitted coefficients
gom_fit$coefficients

         A          B 
0.00010772 0.07534542

Gompertz Model

# Fitted coefficients
gom_fit$coefficients

         A          B 
0.00010772 0.07534542

# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

        A 
0.0021938

Gompertz Model

# Fitted coefficients
gom_fit$coefficients

         A          B 
0.00010772 0.07534542

# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

        A 
0.0021938

# Can simply use the supplied function to calcualte 
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)

$hx
[1] 0.0021938

$par
         A          B 
0.00010772 0.07534542 

$Sx
[1] 0.97269

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)

        A 
0.0021962

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)

        A 
0.0021962

# Compare with life table probability (qx) at age 40:
lt[41,c("age","qx")]

# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1    40 0.00208

Gompertz Model

Gompertz doesn’t fit all age ranges well—particularly young ages.
(This is a well-known fact in demography.)
If we were to focus only on adults, however, this would be a nice way to go.

Gompertz Model

plot(gom_fit)

Gompertz Model

If we were to focus only on adults, however, this would be a nice way to go.

Gompertz Model: 40 Year Old Cohort

ages     <- lt$age[lt$age>=40 & lt$age<100]
deaths   <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]

gom_fit40 <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "gompertz",
                opt.method = "LF2")

Gompertz Model: 40 Year Old Cohort

plot(gom_fit40)

Heligman-Pollard

While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!

Heligman-Pollard

While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP in lieu of gompertz!

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

hp_fit <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "HP",
                opt.method = "LF2")

Heligman-Pollard

plot(hp_fit)

Heligman-Pollard

Mortality rate for a 40 year old:

HP(x = 40, par = hp_fit$coefficients)

$hx
[1] 0.0020409

$par
           A            B            C            D            E           F_ 
 0.000386408  0.021271772  0.107422357  0.001025285  2.871760263 33.201213579 
           G            H 
 0.000022277  1.102611045

x   = 40
mu1 = A^((x + B)^C) + G * H^x)
mu2 = D * exp(-E * (log(x/F_))^2))
eta = ifelse(x == 0, mu1, mu1 + mu2)
hx  = eta/(1 + eta)

Death Hazard at Age 40: Gompertz vs. HP

gompertz(x = 40, par = gom_fit$coefficients)$hx

[1] 0.0021938

HP(x = 40, par = hp_fit$coefficients)$hx

[1] 0.0020409

Alive-Dead (Modeled Mortality)

We’ll Approach This From Two Angles

Discrete Event Simulation: Sample a background mortality time.
Discrete Time Markov: Use parameterized mortality in Alive-Dead model.

1. Sampling Death Dates (DES)

Difficult to sample death times from life tables due to coarseness of data (e.g., five-year age bins).
Solution:
1. Fit a parametric mortality model (e.g., Heligman-Pollard)
2. Construct a survival function from modeled hazards.
3. Sample a death time from the survival function.

1. Sampling Death Rates (DES)

Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)

1. Sampling Death Rates (DES)

Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)

u = 0.25

1. Sampling Death Rates (DES)

We can feed the quantile (0.25) and coefficients into an inverse quantile function.
Straightforward to do for Gompertz model

u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])

[1] 70.467

1. Sampling Death Rates (DES)

Let’s verify this gets us close to what we see in the life table data.

u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])

[1] 70.467

age	mx	Sx
70	0.01827	0.76590
71	0.01995	0.75077
72	0.02162	0.73472

1. Sampling Death Rates (DES)

There is no inverse quantile function for Heligman-Pollard, but we can create one!

1. Sampling Death Rates (DES)

There is no inverse quantile function for Heligman-Pollard, but we can create one!
Approxfun returns a function that provides a value of y for a given value of x.

qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))

1. Sampling Death Rates (DES)

For approxfun():
- x is quantile in CDF for survival, i.e., 1 - Survival¹
- y is the age.
So qHP(0.5) would return median survival age in the modeled data.

1. Sampling Death Rates (DES)

qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)

[1] 71.069

Survival in life table:

age	mx	Sx
70	0.01827	0.76590
71	0.01995	0.75077
72	0.02162	0.73472

1. Sampling Death Rates (DES)

We now have the tools we need to sample death times that closely match underlying mortality data!
For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.

1. Sampling Death Rates (DES)

We now have the tools we need to sample death times that closely match underlying mortality data!
For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.

# Sample death times for 10 individuals
u <- runif(10, min = 0, max = 1)
qHP <- 
  approxfun(y = 0:115,
  x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(u)

 [1] 84.613 77.082 87.010 78.662 86.211 92.977 68.957 89.789 64.050 86.566

2. Modeled Mortality (Markov)

Let’s now incorporate modeled mortality parameters into our Alive-Dead model.
Rather than draw a probability of death from the life table data, we’ll calculate the probability of death at a given age using the model coefficients.

Parameterize the model

params = 
  list(
    t_names = c("lifetable",          # Strategy names
                "heligman-pollard",
                "gompertz"),          
    n_treatments = 3,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 40,                 # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt,                  # Processed HMD life table data (2019, USA)
    gom_fit = gom_fit,                # Fitted Gompertz model object
    hp_fit = hp_fit                   # Fitted Heligman-Pollard model object
  )

Age-Dependent Transition Matrices

fn_mPt <- function(t, params) {
  with(params, {
    h = cycle # Time step
    lapply(t, function(tt){
       # Life table method
        current_age_lt = pmin(initial_age  + (tt)*h - 1,max(lt$age))
        p_death <- lt[lt$age==current_age_lt,"qx"]
        # parametric mortality methods
        current_age = initial_age  + (tt)*h - 1
        r_death_gompertz <- gompertz(current_age,params$gom_fit$coefficients)$hx
        p_death_gompertz <- 1 - exp(-r_death_gompertz)
        r_death_hp <- HP(current_age, params$hp_fit$coefficients)$hx
        p_death_hp <- 1 - exp(-r_death_hp)
        
        array(data = c(1 - p_death, 0, 
                 p_death,1,
                 
                 1 - p_death_hp, 0, 
                 p_death_hp,1,
                 
                1 - p_death_gompertz, 0, 
                 p_death_gompertz,1),
              
              dim = c(n_states, n_states, n_treatments),
              dimnames = list(from = s_names,
                              to = s_names,
                              t_names)) 
    })
  })
}

Construct a Markov Trace

# Add transition matrices into params.
params$mP <- 
  fn_mPt(1:params$n_cycles, params)

# Construct the Markov trace. 
trace <- 
  sim_cohort(params)

Construct a Markov Trace

trace[["lifetable"]][1:10,]

   Alive    Dead
0 100000    0.00
1  99792  207.68
2  99584  416.43
3  99365  634.88
4  99138  861.67
5  98899 1101.19
6  98639 1361.45
7  98358 1641.68
8  98054 1945.62
9  97728 2271.99

trace[["gompertz"]][1:10,]

   Alive    Dead
0 100000    0.00
1  99781  219.14
2  99545  454.89
3  99292  708.46
4  99019  981.15
5  98726 1274.35
6  98410 1589.51
7  98072 1928.21
8  97708 2292.12
9  97317 2682.98

trace[["heligman-pollard"]][1:10,]

   Alive    Dead
0 100000    0.00
1  99796  203.88
2  99584  416.04
3  99363  637.42
4  99131  869.10
5  98888 1112.31
6  98632 1368.43
7  98361 1639.00
8  98074 1925.75
9  97769 2230.56

Calculate Life Expectancy

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr * (1 / params$n_cohort)) %*% payoff_life_exp)

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = lapply(life_exp_cycle,cycle_adj)

Life Expectancy: Alive-Dead Model

total_life_exp

$lifetable
[1] 40.959

$`heligman-pollard`
[1] 41.008

$gompertz
[1] 41.29

Cause-Deleted Mortality

Motivation

Often, background mortality is drawn directly from life table data.
Models also frequently include cause-specific death as a separate tracked event or health state.

Motivation

This approach “double-counts” death because cause-specific death is also reflected in the life table estimates!
OK if cause-specific death is a small contributor to overall death rates.
But what if modeled disease is a frequent cause of death (e.g., cardiovascular disease, cancer)?

Solution

We can construct an alternative life table that “nets out” cause-specific deaths.
Not foolproof: need to make an assumption that we can parse out deaths as indepdendent contributors to overall mortality.

CVD Natural History Model

Next set of slides will walk through the process of constructing a natural history model for cardiovascular disease.
We’ll construct a model with non-CVD background mortality, and CVD mortality.
First step is constructing a cause-deleted life table.
- Needed because at older ages, CVD is a substantial contributor (~50%) of overall mortality.

CVD-Deleted Life Table

Constructing the cause-deleted life table requires two necessary inputs:

Overall mortality by age.
Cause-specific mortality by age (ideally as a % of all deaths at a given age).

These are obtainable for many countries from the Global Burden of Disease and Human Cause of Death Database websites.

CVD-Deleted Life Table

Our estimates of CVD death by age will be drawn from the Global Burden of Disease website.

Screenshot of Global Burden of Disease Data

We extract the percentage of overall deaths that are from CVD by age bin:

ihme_cvd <- 
  tibble::tribble(
        ~age_name,        ~val,
               1, 0.038771524,
               5, 0.038546046,
              10, 0.044403585,
              15, 0.033781126,
              20, 0.035856165,
              25, 0.053077797,
              30, 0.086001439,
              35, 0.130326551,
              40, 0.184310334,
              45,  0.21839762,
              50, 0.243705394,
              55, 0.256334637,
              60,  0.26828001,
              65, 0.272698709,
              70,  0.28529754,
              75, 0.310642009,
               0, 0.016750489,
              80, 0.353518012,
              85, 0.399856716,
              90, 0.447817792,
              95, 0.495305502
        ) %>% 
    mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE))  %>% 
    select(age_ihme,  pct_cvd = val)

Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:

age	age_ihme	pct_cvd	dx	dx_i
0	[0,1)	0.01675	553.221	9
10	[10,15)	0.04440	11.518	1
25	[25,30)	0.05308	103.706	6
50	[50,55)	0.24371	365.721	89
75	[75,80)	0.31064	1997.766	621
98	[95,105)	0.49531	1262.593	625

We next calculate the cause-deleted probability of death (qd) and death rate (md), as well as the cause-specific probability of death (qi) and death rate (mi):

Steps

Cause-specific probability of death (qi): $q \cdot dx_i / dx$
Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.

This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.

Steps

Cause-specific probability of death (qi): $q \cdot dx_i / dx$
Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
- This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.

lt_ <- 
  lt_ %>% 
    # Cause specific probability of death. 
    mutate(qi = q * dx_i / dx) %>% 
    mutate(Rd = (dx - dx_i) / dx) %>% 
    mutate(md = m * Rd) %>% 
    mutate(mi = m - md) %>% 
    mutate(dxd = dx - dx_i)

age	dx	dxd	dx_i	q	qi	md	mi
0	553.221	544.221	9	0.00553	0.00009	0.00546	0.00009
1	38.180	37.180	1	0.00038	0.00001	0.00037	0.00001
2	23.160	22.160	1	0.00023	0.00001	0.00022	0.00001
3	17.689	16.689	1	0.00018	0.00001	0.00017	0.00001
4	14.209	13.209	1	0.00014	0.00001	0.00013	0.00001
5	13.213	12.213	1	0.00013	0.00001	0.00012	0.00001

Modeling CVD-Deleted Mortality

We can also apply the methods covered earlier to the cause-deleted life table data!

ages_     <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_   <- lt_$dxd[lt_$age<100 & lt_$age>=25] 
exposure_  <- lt_$lx[lt_$age<100 & lt_$age>=25]

mort_fit_CVDdeleted <- MortalityLaw(
    x  = ages_,
    Dx  = deaths_,   # vector with death counts
    Ex  = exposure_, # vector containing exposures
    law = "HP2",
    opt.method = "LF2")

Basic CVD Model

Healthy Death based on modeled (Heligman-Pollard) cause-deleted mortality.

Basic CVD Model

CVD Death based on age- and cause-specific death rates from cause-deleted life table.

Basic CVD Model

We also allow for a CVD Death transition based on non-CVD causes.

Age-Dependent Transition Matrices

Construct a Markov Trace

trace <- 
  sim_cohort(params)
trace[["natural_history"]][1:10,]

  Healthy     CVD CVDDeath    Dead
0  100000     0.0   0.0000    0.00
1   90391  9506.2   0.2945  103.00
2   81702 18088.0   1.1404  208.58
3   73847 25833.6   2.4853  316.88
4   66745 32822.6   4.2817  428.11
5   60324 39126.6   6.8544  542.45
6   54519 44809.4  11.5751  660.14
7   49270 49930.7  17.3944  781.42
8   44526 54544.1  23.8195  906.57
9   40236 58697.0  31.3741 1035.87

Good News

Modeled CVD and non-CVD mortality closely matches the cause-deleted life table data!

Name	Model	Code
Infant mortality
Opperman	\(\mu[x] = A/\sqrt{x} - B + C\sqrt{x}\)	opperman
Weibull	\(\mu[x] = 1/\sigma * (x/M)^{(M/\sigma - 1)}\)	weibull
Accident hump
Inverse-Gompertz	\(\mu[x] = [1-\exp((M-x)/\sigma)]/[\exp((M-x)/\sigma)-1]\)	invgompertz
Inverse-Weibull	\(\mu[x] = 1/\sigma * (x/M)^{[-M/\sigma - 1]} / [\exp((x/M)^{(-M/\sigma)}) - 1]\)	invweibull

Name	Model	Code
Adult mortality
Gompertz	\(\mu[x] = A \exp[Bx]\)	gompertz
Gompertz	\(\mu[x] = 1/\sigma \exp[(x-M)/\sigma]\)	gompertz0
Makeham	\(\mu[x] = A \exp[Bx] + C\)	makeham
Makeham	\(\mu[x] = 1/\sigma \exp[(x-M)/\sigma)] + C\)	makeham0
Perks	\(\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]\)	perks
Strehler-Mildvan	\(\mu[x] = K \exp[-V_0 (1-Bx)/D]\)	strehler_mildvan
Adult and/or old-age mortality
Van der Maen	\(\mu[x] = A+Bx+Cx^2 + I/[N-x]\)	vandermaen
Beard	\(\mu[x] = \exp(Bx)/[1+KA \exp(Bx)]\)	beard
Beard-Makeham	\(\mu[x] = A \exp(Bx)/[1 + KA \exp(Bx)]\)	beard_makeham
Gamma-Gompertz	\(\mu[x] = A \exp(Bx)/(1+AG/B[\exp(Bx)-1])\)	ggompertz

Name	Model	Code
Old-age mortality
Van der Maen	\(\mu[x] = A + Bx + I/[N-x]\)	vandermaen2
Quadratic	$= A + Bx + Cx^2%	quadratic
Kannisto	\(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]\)	kannisto
Kannisto-Makeham	\(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]+C\)	kannisto_makeham

Name	Model	Code
Full age range
Thiele	\(\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}\)	thiele
Wittstein	\(q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}\)	wittstein
Heligman-Pollard	\(q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x\)	HP
Heligman-Pollard	\(q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]\)	HP2
Rogers-Planck	\(q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}\)	rogersplanck
Martinelle	\(\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}\)	martinelle
Carriere	\(l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)\)	carriere1
Kostaki	\(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\)	kostaki